#import libraries
library(readr)
library(dplyr)
library(tidyr)
library(tidyverse)
library(ordinal)
library(effectsize)
library(emmeans)

#load wide-format data
data_wide <- read_csv("path", show_col_types = FALSE)
data_wide <- data_wide %>% mutate(id = row_number())
colnames(data_wide)


#transform the wide-format data to the long format
data_long <- data_wide %>%
  pivot_longer(
    cols = c(insurance_high, insurance_low),
    names_to  = c(".value", "inequality"),  # puts 'insurance' & 'risk' into value columns
    names_sep = "_"
  ) %>%
  rename(insurance_option = insurance) %>%   # optional: match your desired column name
  mutate(
    id               = factor(id),
    presentation_order            = factor(presentation_order, levels = c('low_inequality_first', 'high_inequality_first')),
    inequality       = factor(inequality, levels = c("low", "high")),
    insurance_option = factor(insurance_option)  # 3-level categorical outcome
  )
##transform insurance option to be continous for ANOVA
data_long$insurance_option_num <- as.numeric(data_long$insurance_option) -1 


#recode important categorical vars with values
data_wide$gender_coded <- recode_factor(data_wide$gender, 
                                   '0' = 'male', 
                                   '1' = 'female',
                                   '2' = 'non-binary or other',
                                   .default = 'missing')


#get the sample information
##gender
gender_count <- table(data_wide$gender_coded) #male = 55, female = 42, non-binary or other = 3
gender_count 

##age
mean(data_wide$age) #M = 42.54
sd(data_wide$age) #SD = 10.9


#RM ANOVA
##anova
model <- aov(
  insurance_option_num ~ presentation_order * inequality + Error(id/inequality),
  data = data_long
)
summary(model)
##get effect size
effectsize:: eta_squared(model, partial = TRUE)
#######################################################
#Error: id
#                       Df Sum Sq Mean Sq F value Pr(>F)
#presentation_order      1   1.62  1.6200    2.07  0.153
#Residuals               98  76.70  0.7827 
#######################################################

#################################################################
#Error: id:inequality
#                              Df Sum Sq Mean Sq F value  Pr(>F)   
#inequality                     1   2.42   2.420   9.602 0.00254 ** (Eta2:.09)
#presentation_order:inequality  1   2.88   2.880  11.427 0.00104 ** (Eta2:.10)
#Residuals                     98  24.70   0.252                   
###################################################################

##get descriptive
###get raw means
aggregate(insurance_option_num ~ inequality, data_long, mean) #low = .61, high = .83
aggregate(insurance_option_num ~ inequality, data_long, sd) #low = .67, high = .79

aggregate(insurance_option_num ~ inequality * presentation_order, data_long, mean)
aggregate(insurance_option_num ~ inequality * presentation_order, data_long, sd)
#########################################################
#  inequality    presentation_order insurance_option_num
#1        low  low_inequality_first                 0.58 (sd = .61)
#2       high  low_inequality_first                 1.04 (sd = .78)
#3        low high_inequality_first                 0.64 (sd = .72)
#4       high high_inequality_first                 0.62 (sd = .75)
########################################################

###get marginal means & contrasts
emmeans(model, ~ inequality * presentation_order)
####################################################################
#inequality presentation_order    emmean    SE  df lower.CL upper.CL
#low        low_inequality_first    0.58 0.102 155    0.379    0.781
#high       low_inequality_first    1.04 0.102 155    0.839    1.241
#low        high_inequality_first   0.64 0.102 155    0.439    0.841
#high       high_inequality_first   0.62 0.102 155    0.419    0.821
####################################################################

pairs(emmeans(model, ~ inequality * presentation_order))
###########################################################################################
#contrast                                               estimate    SE  df t.ratio p.value
#low low_inequality_first - high low_inequality_first      -0.46 0.100  98  -4.581  0.0001
#low high_inequality_first - high high_inequality_first     0.02 0.100  98   0.199  0.9972
###########################################################################################


#Cumulative Link Mixed Model (CLMM)
##test the proportional odds assumption
prop_odds_test <- clm(insurance_option ~ inequality * presentation_order,
                      data = data_long,
                      link = "logit", threshold = "flexible")
nominal_test(prop_odds_test)
#all ps > .05 -> assumption is met
##fit the regression
fit <- clmm(insurance_option ~ inequality * presentation_order + (1 | id), data = data_long, link = "logit")
summary(fit)
################################################################################################
#Coefficients:
#                                                        Estimate Std. Error z value Pr(>|z|)    
#inequalityhigh                                           1.9887     0.5036   3.949 7.85e-05 ***
#presentation_orderhigh_inequality_first                  0.2370     0.6590   0.360  0.71908    
#inequalityhigh:presentation_orderhigh_inequality_first  -2.1038     0.6942  -3.031  0.00244 ** 
################################################################################################

